Fast cross-staining alignment of gigapixel whole slide images with application to prostate cancer and breast cancer analysis

Joint analysis of multiple protein expressions and tissue morphology patterns is important for disease diagnosis, treatment planning, and drug development, requiring cross-staining alignment of multiple immunohistochemical and histopathological slides. However, cross-staining alignment of enormous gigapixel whole slide images (WSIs) at single cell precision is difficult. Apart from gigantic data dimensions of WSIs, there are large variations on the cell appearance and tissue morphology across different staining together with morphological deformations caused by slide preparation. The goal of this study is to build an image registration framework for cross-staining alignment of gigapixel WSIs of histopathological and immunohistochemical microscopic slides and assess its clinical applicability. To the authors’ best knowledge, this is the first study to perform real time fully automatic cross staining alignment of WSIs with 40× and 20× objective magnification. The proposed WSI registration framework consists of a rapid global image registration module, a real time interactive field of view (FOV) localization model and a real time propagated multi-level image registration module. In this study, the proposed method is evaluated on two kinds of cancer datasets from two hospitals using different digital scanners, including a dual staining breast cancer data set with 43 hematoxylin and eosin (H&E) WSIs and 43 immunohistochemical (IHC) CK(AE1/AE3) WSIs, and a triple staining prostate cancer data set containing 30 H&E WSIs, 30 IHC CK18 WSIs, and 30 IHC HMCK WSIs. In evaluation, the registration performance is measured by not only registration accuracy but also computational time. The results show that the proposed method achieves high accuracy of 0.833 ± 0.0674 for the triple-staining prostate cancer data set and 0.931 ± 0.0455 for the dual-staining breast cancer data set, respectively, and takes only 4.34 s per WSI registration on average. In addition, for 30.23% data, the proposed method takes less than 1 s for WSI registration. In comparison with the benchmark methods, the proposed method demonstrates superior performance in registration accuracy and computational time, which has great potentials for assisting medical doctors to identify cancerous tissues and determine the cancer stage in clinical practice.

www.nature.com/scientificreports/ drug development, IHC has been applied to test drug efficacy by monitoring certain protein expressions and analyzing the active levels of disease targets 5,6 . In clinical diagnosis, with the knowledge of tumor marker and IHC tools, medical professionals and scientists are able to diagnose tumors as benign or malignant, determine the cancer stage and identify the cell type and origin of a metastasis in order to find the site of the primary tumor. However, to align, compare, and quantify protein expressions of tissues of interest from multiple whole slide images (WSIs) with different stainings is difficult due to the enormous dimensions of gigapixel images to analyze, nonrigid distortions and deformations introduced during slide preparation and dissimilar cell appearance and tissue morphology across stainings. Figure 1a shows challenges induced from individual data preparation steps for cross staining WSI analysis in cancer diagnosis. Serial sections are firstly cut from a patient donor block and mounted onto glass slides, causing nonrigid deformations and distortions across slides. Secondly, slides are processed with various kinds of stainings. As each staining highlights specific substances in tissues, the image appearances of cells and tissue layouts across various stainings become distinctively dissimilar. Thirdly, slides are digitized into enormous gigabyte-sized WSIs at single-cell microscopic resolution, including a H&E WSI and multiple IHC WSIs to observe tissue morphology and various protein expression patterns together. Typically, WSIs are stored in a gigabyte-sized multi-resolution pyramid data structure from low to high magnification with a general size of around 100,000 × 200,000 pixels at the highest resolution level. In addition to the enormous dimensions of data to process, there are large variations on stain colors, cell appearance, and tissue morphology across WSIs.
Due to the advancement of computing power, many computer algorithms [7][8][9][10] have been developed for automatic alignment of biological microscopic images. The bUnwarpJ tool by Arganda-Carreras et al. 7 is developed for biological section alignment, allowing bi-directional image registration. Saalfeld et al. presented a least square based registration method using linear feature correspondences 8 and an elastic registration approach based on non-linear block correspondences 9 . Wang et al. built a CwR method 10 , which contains a 3D alignment and validation model utilizing the b-spline deformation field. In this study, the above mentioned approaches 7-10 are adopted as benchmark methods.
From the literature review, previously many studies [11][12][13][14][15][16][17][18][19][20] performed image alignment by directly applying image transformation to entire images. Later, several approaches [21][22][23][24][25][26] have been developed with a more efficient scheme such as hierarchical frameworks to deal with large microscopic datasets. Due to the advancement of imaging technology, image alignment is required in application to enormous gigabyte-sized images such as WSIs. The above-mentioned methods are however too memory costly and slow to be applied to WSI alignment. For dealing with the enormous size of WSIs, several tile-based approaches have been built. Roberts et al. 27 presented a multi-resolution-block-matching-based nonrigid registration framework. Song et al. 28 and Lotz et al. 29 extended Roberts's approach 27 for tissue reconstruction of histological sections with different stains. Song et al. developed an unsupervised content classification method that produces multichannel probability images from images of different stains to increase the similarity of the two images and integrated the classification method into the multi-resolution-block-matching-based framework. Lotz et al. 29 first computes nonlinear registration on low-resolution images to correct global large-scale deformations occurring in tissues, and the results of this nonlinear registration is then used as an initial guess for a patch-wise registration. In this study, the three tilebased approaches for WSI alignment [27][28][29] are adopted as benchmark methods for run time analysis.
Prostate cancer is the second leading cause of cancer related death in men in the United States 30 . In 2018, prostate cancer has caused 358,989 deaths worldwide [31][32][33] . In diagnosis of invasive prostate carcinoma, it is vital to check if the basal cell layer is absent or not. Hence, a complete absence of staining in basal cell-associated markers is supportive of a malignant interpretation. The diagnostic and prognostic value of CK18 in prostate cancer and human tumor has been proved and described in [34][35][36] . In 2012, Weng et al. 34 pointed out that CK18 is expressed in a variety of adult epithelial organs, plays a role in carcinogenesis and in combination with other markers as an important marker of therapeutic efficacy. In 2016, Yin et al. 35 determined that the IHC CK18 expression in PCa tissues is inversely correlated with tumor grade in a statistically significant fashion ( p = 0.028 ) where tissue samples with higher tumor grades (Gleason score ≥ 7 ) show gradually decreased CK18 intensity. In 2021, Menz et al. 36 show that reduced CK18 immunostaining is linked to high UICC stage ( p = 0.0010 ), high Thoenes grade ( p = 0.0086 ) and advanced tumor stage ( p < 0.0001 ) based on 11952 tumor samples. Using IHC, High-Molecular-Weight Cytokeratin (HMCK) is a cytoplasmic marker that highlights intermediate cytokeratin (CK) filaments in glandular basal cells and is specific for basal cells in the prostate.
In clinical practice, manual cross-staining analysis of whole slides is conducted; in order to diagnose a sample as benign or malignant, medical doctors have to manually identify cancerous tissues first and then determine the tumour stage and grade of a patient with prostate cancer. As only the tissues with positive CK18 detected and within the gland with positive HMCK detected are determined as cancerous tissues, in order to identify cancerous tissues and determine the tumor stage and grade, medical doctors have to manually align three gigabyte-size microscopic WSIs of each patient and then subjectively quantify cancerous tissues of WSIs, to perform simultaneous analysis of an IHC HMCK slide, an IHC CK18 slide and a H&E slide. This is however difficult and poorly reproducible considering the enormous dimensions of gigapixel images to analyze, nonrigid distortions and deformations introduced during slide preparation and dissimilar cell appearance and tissue morphology across different staining. Automatic cross-staining alignment of WSIs enables medical experts and scientists to easily identify the cancerous tissues, quantify the amount of cancerous tissues and determine the stage and grade of a tumour. Figure 1b shows the automatic triple-staining result of prostate cancer sample by the proposed method.
Breast cancer is the most prevalent cancer in women and a major cause of mortality 37,38 . In prognosis of patients with breast cancer, the five-year survival rate is 90%, and ten-year is 83%. However, the survival rate is significantly worsened when breast cancer metastasizes 38,39 . In the case of localized breast cancer, the five-year survival rate is 99%, which drops to 85% in the case of regional (lymph node) metastases. Furthermore, the fiveyear survival rate is 26% in the case of distant metastases. Therefore, it is important to identify the metastases to www.nature.com/scientificreports/ provide adequate treatment and to improve the survival rate. The prognosis of breast cancer patients is mainly determined by the extent of the cancer evaluated using the TNM staging system that assesses the size of the tumour (T), the status of axillary lymph nodes (N) and the status of metastasis (M). The status of axillary lymph  www.nature.com/scientificreports/ nodes could be estimated by evaluating the sentinel lymph node, which is the first lymph node to receive the afferent lymphatic drainage from the primary cancer and is thus the first node to be involved by metastasis. Therefore, patients whose sentinel lymph node is negative for breast cancer metastasis could be spared for more extensive axillary lymph node dissection. Alternatively, axillary lymph node dissection may be performed on breast cancer patients, particularly for patients with high suspicion of axillary lymph nodes metastasis. Metastatic foci to lymph nodes have been classified into macrometastases (metastatic size greater than 2.0 mm), micrometastases (metastatic size greater than 0.2 mm, but none greater than 2.0 mm), and isolated tumour cells (ITCs, metastatic size no greater than 0.2 mm). It is difficult for human to identify ITCs due to their tiny size, particularly in cases with many axillary lymph nodes, leading to under-treatment of the patients. Furthermore, the routine examination of H&E slides of lymph nodes, particularly for axillary lymph node dissection, which may contain many lymph nodes, is time consuming and pathologists are at the risk of missing tiny metastatic foci. The MIRROR (Micrometastases and Isolated tumour cells: Relevant, Robust or Rubbish) trial emphasized the need for additional therapy in cases of invasive breast carcinoma with ITCs or micrometastases [40][41][42][43] , supporting the identification of ITCs or micrometastases is urgently needed in the management of breast cancer patients. IHC staining for cytokeratin on lymph node section could be used to detect ITCs. In this study, we illustrate automatic cross-staining alignment of WSIs of H&E and IHC CK(AE1/AE3) slides help screen all the lymph nodes and identify tiny metastatic foci to improve the accuracy of pathological assessment of lymph node as shown in Fig. 1c.
In this paper, we present a real time interactive fully automatic hierarchical registration framework that is able to perform cross-staining alignment of gigapixel WSIs of histopathological and immunohistochemical microscopic slides in real time and also overcome the major speed bottleneck for assisting medical doctors to identify cancerous tissues and determine the stage and grade of a tumor. The proposed hierarchical registration framework is described in detail in "Methods" section . The effectiveness and robustness of the proposed fully automatic hierarchical registration framework is validated using two datasets, including a dual staining breast cancer dataset with 43 H&E slides and 43 IHC CK(AE1/AE3) slides and a triple staining prostate cancer dataset with 30 H&E slides, 30 IHC CK18 slides, and 30 IHC HMCK slides. Figure 1b,c presents the automatic WSI alignment result of triple-staining prostate cancer samples and dual-staining breast cancer samples by the proposed method, respectively. Figure 2 illustrates the flowchart for the proposed method in cross staining WSI analysis in clinical usage. In addition, an online web-based system of the proposed method has been created for live demonstration of real time cross staining alignment of whole slide images. Please see the supplementary video at https:// www. youtu be. com/ watch?v= 0Uc6-s_ ClIg& ab_ chann el= ProfC hing-WeiWa ng.  The slides were generated under routine sample preparation procedures. Firstly, tissues are fixed in 10% buffered formalin and embedded in paraffin. Next, serial sections were cut (4 μm) using Leica RM2155 and then stained with H&E and associated IHC stains. In digitization, the prostate cancer slides were scanned using Leica AT Turbo scanner (Leica. Wetzlar Germany) at 40× objective magnification and the breast cancer slides were scanned using 3DHISTECH Pannoramic SCAN II scanner (3DHISTECH Kft. Budapest Hungary) at 20× objective magnification. We summarized the information of the two experimental data sets in Table 1.

Results.
In evaluation, three experiments were performed, including an experiment on cross-staining alignment on triple staining prostate cancer slides (see "Quantitative evaluation on the triple-staining prostate cancer dataset" section ), an experiment on cross-staining alignment on dual staining breast cancer slides (see "Quantitative evaluation on the dual-staining breast cancer dataset" section ) and an experiment on run time analysis (see "Run time analysis" section ). The first evaluation was conducted on 30 sets of triple-staining prostate cancer samples, including H&E slides and IHC slides with two kinds of antigens, CK18 and HMCK, and we compare the proposed method with four image registration techniques for biological microscopic image alignment, including bUnwarpJ 7 , LeastSquares 8 , Elastic 9 and CwR 10 . The second evaluation was performed on 43 pairs of dual-staining breast cancer slides, including 43 H&E slides and 43 IHC slides with CK(AE1/AE3) antigen, and we further compare the proposed method with the best benchmark approach performed in the first experiment, i.e bUnwarpJ 7 . The third experiment was performed to produce run time analysis on the 43 pairs of dual-staining breast cancer slides. We compare the computational time in WSI registration of the proposed method and the best benchmark approach performed in the first experiment, i.e bUnwarpJ 7 , Roberts et al. 27 , Song et al. 28 and Lotz et al. 29 . All statistical analysis was performed using SPSS software 44 , and all experiments were performed on a workstation with an Intel Xeon Gold 6134 CPU processor and 64GB RAM. Regarding the evaluation method to measure the registration accuracy, as pointed out in the previous studies 45,46 , the traditional evaluation approach based on sum of squared differences of image intensities between the target and the transformed source could be inaccurate in biological image applications as the pixel intensity in the target and the intensity of the accurately registered pixel in the transformed source are generally different due to stain variation. As a result, the quantitative evaluation approach of the previous works 45,46 is adopted where five corresponding landmarks between target images and associated transformed source images by each registration method were firstly manually marked by experienced pathologists, and an automatic matching system is applied to compare the coordinates of the corresponding manual annotations, producing registration accuracy scores for each image pair based on the matching successful rate over the corresponding annotations (within a five pixel distance). The final registration accuracy score of each registration method is computed using the averaged accuracy over all image pairs as shown in Fig. 3.
Quantitative evaluation on the triple-staining prostate cancer dataset. Figure 4a presents the evaluation results of registration accuracy on triple-staining prostate cancer tissue images by the proposed method and four benchmark approaches, including bUnwarpJ 7 , Leastsquares 8 , Elastic 9 and CwR 10 . In addition, four transformation models, including affine, rigid, similarity and translation, are applied in testing of the benchmark approaches except for the CwR 10 . As shown in Fig. 4, the proposed method consistently performs well and significantly outperforms benchmark methods. Detailed quantitative results could be found in Table 2. In addition, we further conduct the statistical analysis, and the result shows that the proposed method is significantly better than the benchmark methods ( p < 0.01 ) for triple-staining prostate cancer tissue image registration. Figure 3 shows some sample results of registration outputs of the proposed method and benchmark approaches.
Quantitative evaluation on the dual-staining breast cancer dataset. Figure 4b compares the evaluation results on dual-staining breast cancer tissue images by the proposed method and the best benchmark approach performed in the first experiment, i.e, bUnwarpJ 7 . The proposed method achieves high averaged registration accuracy (93.49%) while the bUnwarpJ approach 7 obtains averaged accuracy only 25.58%. Detailed quantitative results www.nature.com/scientificreports/ www.nature.com/scientificreports/ could be found in Table 3. Using SPSS software 44 , the result shows that the proposed method significantly performs better than bUnwarpJ 7 ( p < 0.001).

Run time analysis.
For run time analysis, the breast cancer dataset is used, and the computational time of WSI registration of the proposed method is compared with the best benchmark approach performed in the first experiment, i.e. bUnwarpJ approach 7 , and three additional benchmark approaches, including Roberts et al. 27 , Song et al. 28 and Lotz et al. 29 . In testing of WSI registration, as the system of bUnwarpJ approach 7 tends to fail due to out of memory, we therefore perform a regression analysis to estimate the computing time for WSI registration using least squares. Detailed quantitative results could be found in     Table 3. Results on the breast cancer dataset in comparison of the proposed method and the best performed benchmark approach in the first experiment. *The proposed method is significantly better than the benchmark approach ( p < 0.001) www.nature.com/scientificreports/ those three different slides, which is time-consuming, subject and difficult to quantify tumors. Automatic crossstaining biological image alignment helps medical doctor to identify cancerous tissues, quantify the amount of cancerous tissues and determine the cancer stage. Figure 1b shows triple-staining alignment on prostate cancer tissue sample results performed by the proposed method. In prognosis of breast cancer patients, metastatic foci to lymph nodes have been classified into macrometastases (metastatic size greater than 2.0 mm), micrometastases (metastatic size greater than 0.2 mm, but none greater than 2.0 mm), and isolated tumour cells (ITCs, metastatic size no greater than 0.2 mm). Currently, the status of lymph nodes is routinely determined by the examination of H&E slides. Immunohistochemical staining for cytokeratin on lymph node section may be applied to rule out ITCs, which are mostly detected by IHC due to their tiny size. Despite the clinical significance of ITCs or micrometastases is debated, the results of the MIRROR trial emphasized the need for additional therapy in cases of invasive breast carcinoma with ITCs or micrometastases [40][41][42][43] , supporting that the identification of ITCs or micrometastases is currently clinically needed. However, because IHC for cytokeratin is not routinely performed on every axillary lymph node of breast cancer patients, it is possible that ITCs may run undetected, particularly in cases with many axillary lymph nodes, and result in under-treatment of the patients. Furthermore, the routine examination of H&E slides of lymph nodes, particularly for axillary lymph node dissection which may include many lymph nodes, is time-consuming and pathologists may potentially run the risk of missing tiny metastatic foci. In this study, automatic dual-staining WSI alignment of H&E and IHC CK(AE1/AE3) slides is performed, which help medical doctors to screen all the lymph nodes sections and identify potential tiny metastatic foci as shown in Fig. 1b. Due to the budge limitation, immunohistochemistry for CK can be performed on two sentinel lymph nodes (LN), usually sentinel LN1 and LN2. However, in clinical practice, we are very frequently encountering more than two sentinel LNs plus axillary lymph node dissections, which may include many lymph nodes. This method can help identifying potential micrometastasis and ITCs that may go unnoticed by the pathologists in lymph nodes sections without immunohistochemistry for CK.

95% confidence interval for mean
In the experiments of triple-staining alignment on 90 prostate cancer tissue images, the proposed method achieves the triple-staining registration accuracy 0.833±0.0674, while four benchmark approaches 7-10 perform poor, obtaining the averaged accuracy no higher than 0.60. In statistical analysis, the proposed method is significantly better than all benchmark approaches 7-10 ( p < 0.01 ) for cross-staining microscopic image registration. In alignment of breast cancer tissue images evaluation, we compare the proposed method and the best benchmark approach performed in the triple-staining alignment of prostate cancer tissue images, i.e bUnwarpJ 7 on 86 breast cancer tissue images. In evaluation, the proposed method achieves high image registration accuracy 0.931±0.0455. In comparison, the benchmark method -bUnwarpJ approach only obtains averaged accuracy 0.255. In statistical analysis, the proposed method significantly outperforms bUnwarpJ approach ( p < 0.001 ) for dual-staining alignment of breast cancer tissue images.
For run time analysis, the proposed method takes only 4.34 s per WSI registration on average. In comparison, the benchmark approaches 7,27-29 spend more than 70 s per WSI registration on average. In addition, for 30.23% data, i.e. 13 out of 43 WSI pairs of breast cancer data, the proposed method takes less than 1 s for WSI registration. In this study, we present a real-time, fully automatic and robust cross-staining registration system for aligning multiple IHC slides and histopathological H&E slides to assist assessment of tissue morphology and various protein expressions together, using two different digital scanners, i.e. Leica and 3DHISTECH, on two different cancer samples, i.e. breast cancer and prostate cancer, and from two different hospitals. The proposed method is not limited for cross-staining analysis but could also be used for single-staining serial section comparison.

Methods
In this study, we present a real-time and fully automatic coarse to fine propagated tile-based registration framework for cross-staining WSI alignment. The proposed framework consists of two main parts, including a rapid global registration module ("Rapid global image registration" section ) and a propagated real time multi-level image registration module ("Propagated multi-level image registration" section). Fig. 5 presents the flowchart of the proposed method. For the rapid global registration module as illustrated in Fig. 5a., firstly, cytoplasm features of low-level images are extracted using color deconvolution ("Cytoplasm feature extraction model" section); secondly, corresponding landmarks between the target and the source image are detected ("Corresponding landmark detection" section) and then used to compute the global transformation parameters ("Computation of global transformation parameters" section); thirdly, the low-level source image are aligned to the target by using the global transformation parameters ("Global image transformation" section). For the second part of the proposed method, i.e. the real time propagated multi-level image registration module as illustrated in Fig. 5b), firstly, the tile set of the source WSI, corresponding to the Field of View (FOV) of the target WSI, are extracted using the transformation parameters; secondly, the source FOV is aligned to the target FOV by the proposed real time patch-wise registration ("Propagated multi-level image registration" section).
All methods were carried out in accordance with relevant guidelines and regulations. The experimental protocols were approved by the research ethics committees of the National Taiwan University Hospital (with ethical approval number: NTUH-REC 201810082RINB) and the Tri-Service General Hospital, Taipei, Taiwan ((with ethical approval number: TSGHIRB1-107-05-171). Informed patient consent forms were formally waived by the research ethics committees of the National Taiwan University Hospital and Tri-Service General Hospital, and the data were de-identified and used for a retrospective study without impacting patient care.
Rapid global image registration. Let { J ξ } be a pair of digital WSI, containing a target WSI J t and a source WSI J s for alignment where ξ = t represents the target, and ξ = s represents the source, respectively. We formulate a set of digital WSIs { J ξ } into multi-level pyramid tile-based data structure { I To obtain the global transformation parameters and low-level image registration result, low-level image registration process comprises four parts as shown in Fig. 5a. (1) Using colour deconvolution technique with an orthonormal transformation matrix, stain separation is performed on the low-level target and source tile-images to extract the cytoplasm features of the target and source. (2) Corresponding landmark detection is applied to the Then, colour c is transformed to the new unit vector � c = u.� u + d. � d + n.� n + � p . By setting u = 0 , the undesired component is removed, generating new colour � c ′ = d. � d + n.� n + � p . In three channels image, a colour system is described as a matrix form µ with every row representing a specific stain and every column representing the Optical Density (OD) as detected by c 1 , c 2 and c 3 for each stain, respectively.
In this study, a normalized OD matrix, µ , to describe the colour system for orthonormal transformation is defined as follows: Let G = [g 0 , g 1 , g 2 ] be denoted as a 3 × 1 vector for the amount of stains at the specific pixel, where g 0 , g 1 and g 2 represent the first, second and third channel, respectively. The OD levels at individual pixels could be formulated as O = G µ . Therefore, multiplication of the OD image with the inverse of the OD matrix results in orthogonal representation of the stains forming the image G = µ −1 O . Then, the cytoplasm features, g 1 0 and g 2 0 , of the target image and the source image are extracted for further corresponding landmark detection.
Let P = { p k } K k=0 and Q = { q k } K k=0 be two sets of pairwise correspondences between g 1 0 and g 2 0 in Ŵ , respectively, where α k and β k are not negative coefficients.
Two sets of pairwise correspondences P and Q are obtained by the following steps: (1) key point sets E 1 and E 2 are detected by applying Difference of Gaussian (DoG) detector 48 , (2) the corresponding landmarks P and Q are selected as geometric consensus between keypoints E 1 and E 2 by applying Random Sample Consensus (RANSAC) 49 . The corresponding feature detection P and Q are then used to compute a set of global transformation parameters. www.nature.com/scientificreports/ Computation of global transformation parameters. By aligning two sets of pairwise correspondences, P = { p k } K k=0 and Q = { q k } K k=0 , the task is to compute a set of global transformation parameters rapidly at low resolution level, including scaling factor S, the rotation matrix R and the translation vector T. Firstly, to computes the scaling factor S, the centroids of each set are calculated, then both of the centroids are translated to its origin (the centred vectors). Let p = K k=0 p k K and q = K k=0 q k K be the centroids of each set, and p ′ k = p k − p and q ′ k = q k − q be the centred vectors, the scaling factor S can be computed as follows: where σ represents the variance. Next, d × d covariance matrix C is computed: , and tr is matrix transpose operator. Then using Jacobi Singular Value Decomposition (SVD) algorithm 50 δ , the SVD matrix X is computed: X = δ(C) . To obtain the rotation matrix R, the SVD matrix X is decomposed, UZV = X , where U, V are rotation matrices and Z is a diagonal matrix. The rotation matrix R can be computed by the following expression: After computing scaling factor S and rotation matrix R, the translation vector T can be computed by: S, R and T are then used for aligning low-level source image I s low to the target I t low , real-time interactive FOV localization and real-time multi-level image registration processes.
Global image transformation. After obtaining a set of transformation parameters from the previous section, next step is to align the low-level source image I s low,(x,y,w,h) to the target image and to produce a low-level registered-source I s ′ low,(x ′ ,y ′ ,w,h) , where (x, y) is the coordinate of I s low,(x,y,w,h) and (x ′ , y ′ ) is the coordinate of I s ′ low,(x ′ ,y ′ ,w,h) . The mapping relationship is formulated as follows.
Propagated multi-level image registration. Real-time interactive FOV localization module is devised to locate and fetch the associated tile set of the source WSI corresponding to the FOV of the target WSI as shown in Fig. 5b. Let F t = I t l,(u t ,v t ,w t ,h t ) be the FOV of the target WSI, where (u t , v t ) is the global coordinate of F t ; n 1 × n 2 is the number of tiles containing F t . Firstly, the global coordinate (u s , v s ) of the top left tile of the FOV in the source WSI corresponding to the target FOV F t is computed, where △ = l − l low is used for level calibration.
Secondly, the associated tile set of the source WSI corresponding to F t is fetched: G : { g l,i,j } { i=a,...,a+n 1 ,j=b,...,b+n 2 } . The associated tile set is then used to compute registration outputs. Next, transformation is applied to each enlarged area of every tile g l,i,j in the fetched tile set, producing a transformed enlarged area. Afterwards, the top left tile of each transformed enlarged area is selected and then integrated into an image as the registration output, i.e. the corresponding FOV of the source WSI (see Fig. 5b). Detailed registration processes of the proposed method are described as follows.

Supplementary Video
An online web-based system of the proposed method has been created for live demonstration of real time cross staining alignment of whole slide images. Please see the supplementary video at https:// www. youtu be. com/ watch?v= 0Uc6-s_ ClIg& ab_ chann el= ProfC hing-WeiWa ng.